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We reexamine the production of gravitational waves by bubble collisions during 
a first-order phase transition. The spectrum of the gravitational radiation is deter- 
mined by numerical simulations using the " envelope approximation" . We find that 
the spectrum rises as f^'^ for small frequencies and decreases as f~^'^ for high fre- 
quencies. Thus, the fall-off at high frequencies is significantly slower than previously 
stated in the literature. This result has direct impact on detection prospects for 
gravity waves originating from a strong first-order electroweak phase transition at 
space-based interferometers, such as LISA or BBO. In addition, we observe a slight 
dependence of the peak frequency on the bubble wall velocity. 



I. INTRODUCTION 

Colliding bubbles in a first-order phase transition constitute one possible source of 
stochastic gravitational wave (GW) radiation If the electroweak phase transition 

is strongly first-order, for instance, the kinetic energy stored in the Higgs field and the bulk 
motion of the plasma is partially released into gravity waves. This happens mostly at the end 
of the phase transition, when collisions break the spherical symmetry of the individual Higgs 
field bubbles. This possibility was systematically analyzed in a series of papers [s], [4, ^,1^. 
The first simulation [sl consisted hereby of the full scalar field dynamics of two bubbles in 
vacuum, where the essential observation was made that the emitted radiation depends only 
on the gross features of the problem, namely the kinetic energy stored in the uncoUided 
bubble regions. This observation is the basis of the so-called envelope approximation that 
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opened up the possibility of simulating phase transitions with a large number of bubbles. 



This was subsequently exploited in refs. 



in ref. jsj. 

In the case of only two colliding bubbles, the spectrum decreases as /~^'^ for high fre- 
quencies. But this result might be special to the case of two bubbles, where the collision 



5| and further refined for a thermal environment 



The 



never finishes, which makes the introduction of a time cutoff function mandatory 
situation is different if a realistic phase transition with a large number of bubbles is sim- 
ulated. There were hints in refs. that the spectrum of multi-bubble simulations 
might be more flat than in the two bubble case, but the numerical accuracy prohibited a 
conclusive statement. As a result, the frequency fall-off of the two bubble case is still being 
used in the present day literature 

The aim of the present work is to reexamine the generated spectrum of GWs by simulating 
a phase transition with a large number of bubbles, making use of the aforementioned envelope 
approximation. Compared to ref. , the numerical accuracy will be considerably improved, 
and a larger portion of the spectrum will be determined to allow for a careful analysis of the 
high frequency behavior. 

II. DETERMINATION OF THE GW SPECTRUM 

The fundamental quantity that enters the determination of the gravitational radiation are 
the spatial components of the stress-energy tensor Tjj(x, t). For a thermal phase transition 
it consists of the scalar field part and the plasma contribution. The total energy radiated 



into a direction k is then given by [1( 

2Gu;^A^j^imik)T*Ak, a;)Tz^(k, to), (1) 



dEcw 



where Tjj(k, a;) denotes the stress-energy tensor in Fourier space 

T,,(k,w) = ^J dte''^' j d'xe-''^^-^T,,i^,t), (2) 
and A is the projection tensor for the transverse-traceless part 

A,,,,^(k) = 6uS,^ - 2k,k^Su + -k,k,k,k„ - -6,,6i^ + -%k,k™ + -5,„k,k,. (3) 
Eq. ([1]) is derived in the wave zone approximation that is well justified in the present case. 



The key idea of the envelope approximation js] is that the GW production does not 
depend on the details of the evolution of the scalar field in the region of the collided bubbles, 
but rather on the gross features of the problem, namely the shape of the uncoUided bubble 
walls. Denoting the times and positions of the nucleated bubbles by t„ and x„, this turns 
eq. (ED into 



Tij(k,uj) 



1 

2^ 



dt e 



(4) 



where S'„ denotes the uncollided region of the nth bubble and Tjj „ its stress-energy tensor. 
If the phase transition proceeds by detonation, the stress energy is concentrated in a shell 
of bulk motion that is thin compared to the bubble radius, thus motivating the thin-wall 



approximation 



An 



(5) 



with Rn{t) the size of the nth bubble and p^ac the difference in energy density between the 
true and the false vacuum (we use this terminology even though in the presence of a thermal 
bath the latent heat is the relevant quantity). 

The efficiency factor k determines how much of the vacuum energy is transformed into 
kinetic energy of the bulk fluid instead of reheating the plasma inside the bubble. This 
coefficient can be determined as a function of the parameter a that is defined as the ratio 
between the vacuum energy and the total energy stored in radiation 
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(6) 



The last ingredient from the phase transition is the velocity of the bubble wall Vb- For very 
strong phase transitions it is given by jll | 



Vb 



\/V3 + 



'a" 



2a/3 



1 + a 



(7) 



In fact, eq. ([7]) is not true in general. In phase transitions, there exists a larger class of 
detonation solutions, as discussed in ref. 121]. Therefore eq. (JTj) gives only a lower bound on 
the wall velocity and we will treat Vb as a free parameter in our analysis. 
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Following the above considerations, eq. (jlj) can be rewritten as 

Tij(k,Uj) = Kp,,acVlCij(k,Uj), (8) 

A„,y(k,cu) = f rffi e-^^^''(*-*")'^-%%. (10) 



An efficient way of evaluating these integrals numerically is presented in the Appendix. 

Before we discuss the simulation of the phase transition, we comment on the so-called 
quadrupole approximation that is given by the limit k ■ x — or 

C,,{u) = f dte''^\t-tnfAn,j{0,t). (11) 

This approximation becomes exact in the limit of small frequencies. For large frequencies, it 
turns out that this substitution in the case of two bubbles overestimates the resulting GW 
production by more than one order of magnitude; additionally, the peak frequency is shifted 
to larger values what was ffist noticed in the simulation of the scalar field evolution of two 
single bubbles in ref. [3|. This result is rather surprising, since e.g. in electrodynamics with 
a localized source, the quadrupole approximation generally underestimates the full result. 
However, in the present case the assumption of localized sources is not applicable such that 
the quadrupole approximation can overestimate the full result as detailed in ref. 

To model the phase transition, we assume that the nucleation probability per volume and 
time is given by 

P = Poexp{P{t-to)). (12) 

The parameter is approximately the duration of the phase transition. Since (3 is the 
sole dimensionful parameter of the phase transition, the fraction of energy liberated into 
gravitational wave radiation per frequency octave is 6j 



where we used the definition of the Hubble constant 

2 _ ^TlGptot _ 8nG{pyac + prad) 

~ 3 ~ 3 ' ^ ' 



and defined the dimensionless function A as 



^{u/l3,v,) = -^^Sr I dkA,,,i^C:^Cim. (15) 
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If the integration is performed over a large volume V of the Universe (in terms of (3 ^) 
with N bubbles, one obtains in the quadrupole approximation 



The remainder of the present work is devoted to discussing numerical results for the 
function A(x, Vb). 



To calculate the function A{x,Vb), we proceed as follows. First, we choose the main 
parameters of the simulation: The radius of the part of the Universe under consideration 
Lu, and the number of directions to integrate over, Nk- Since the GW production is isotropic, 
the directional variation indicates the statistical uncertainty of the simulation. 

Using this information, we simulate a scenario with random nucleation times and locations 
for the bubbles following the nucleation probability in eq. f|T2|) . As in ref. our volume is 
spherical and the bubbles are cut off by the boundary which is equivalent to having a mirror 
symmetric configuration beyond the volume. 

The result for a small and a large Universe and a few directions, Nk = 32, is shown 
in Fig. [H For large wall velocities, the size of the simulated part of the Universe mostly 
influences the total amplitude of the GW radiation. However, for smaller wall velocities, 
the produced radiation at high frequencies depends crucially on the size of the Universe. 
For small velocities the high frequency part of the spectrum approaches the quadrupole 
approximation, since in this case k • x k ■ Xn. This is not identical to the quadrupole 
approximation (only the phase in eq. ( ITOi) vanishes), but it turns out that the additional 
phase in eq. iQ from the bubble position solely leads to a slight suppression of the low 
frequency part of the spectrum in our numerical results. 

The enhancement of the high-frequency part in the quadrupole approximation can be 
explained in the following way. After the nucleation of the first bubble, the phase transition 
ends when the whole space is in the new phase, which typically takes a time r ~ 5/(3. 




(16) 



and hence we recover the well known result 6| 




(17) 



III. NUMERICAL RESULTS 
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FIG. 1: The left panel shows the fraction of gravitational radiation for a small simulation, Ljj = 
3vb/(3, 7 bubbles. For small wall velocities, significant boundary effects are visible. The right 
simulation is relatively large, Ljj = 7vb/P, 109 bubbles. Both simulations have been integrated 
over Nk = 32 directions. Velocities decrease from top to bottom as Vb = {1, 0.1, 0.01}. 



However, if the radius of the simulated volume Lu is chosen too small, the simulation 
ends when the first nucleated bubble covers the whole spherical volume. This leads to a 
suppression of the total GW density and also to a rather steep edge in the time profile 
that amplifies the high frequency part of the spectrum. This enhancement effect is non- 
physical and just arises because of the finite size of our test volume and the presence of the 
boundary. It disappears for larger volumes, when the phase transition is finished before the 
first nucleated bubble covers the whole spherical volume and boundary effects are small. For 
large bubble velocities, when the quadrupole approximation does not apply, the enhancement 
also disappears, since effects from large bubbles are in general reduced due to the fast 
oscillation of the integrand in eq. ffTOl) . This motivates our final choice Lu = 7vb/P in order 
to suppress these boundary effects. 

The error bars in Fig. [T] result from the integration over the Nk = 32 different directions. 
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If the number of directions would be further increased, most of the displayed oscillatory 
features in the spectrum would persist, since they result from the finiteness of the size of the 
simulation Lu and depend on the specific nucleation positions and times of the randomly 
generated scenario. In addition, the amplitude of the low frequency part of the spectrum and 
the total emitted radiation depend sensitively on the times of the first bubble nucleation and 
first collision, which are however quite similar in the different scenarios. To decrease these 
features, one has to further increase the size of the Universe and/or average over several 
scenarios. We will do the latter in the following. 

Our final results, shown in Fig. [21 are obtained for the parameters 

Lu = 7vb/P, Nk = 32, (18) 

and averaged over eight scenarios. The given error bars result from this averaging procedure. 
For very small frequencies, one expects flaw* to scale as u^, since Cij{u) becomes constant 
in this limit. This is reflected by the numerical results. For large frequencies, we find 



that ^GW* scales as uj ^, unlike the two-bubbles case that was discussed in ref. js]. We 
parametrize the spectrum as 

0/* ~r CtJ* 

with the peak frequency /* = u;/27r and the peak amplitude CIqw* and both parameters 
are functions of the wall velocity. The two exponents obtained by fitting lie in the range 
a G [2.66,2.82] and b G [0.90,1.19]. The most interesting case from a phenomenological 
point of view is given by large wall velocities, Vb ~ 1, in which case the fit yields a ^ 2.8 and 
b ^ 1.0. Note that the fit is optimized for a frequency range close to the peak frequency, 
and does not correctly reproduce the asymptotic low frequency behavior. 

The main result of our analysis is that the spectrum scales for high frequencies close to 
uj~^ in contrast to the two bubble case presented in ref. 3l]. This feature is already present 
for a rather small volume (just a few bubbles) and high wall velocities {vb — 1) as depicted 



in the left panel of Fig. [TJ The reason for this qualitative difference between t 



le results is 



hence probably the time-dependent cutoff function that has been used in ref. jsl in order to 
terminate the phase transition. 

The spectrum that would be observed today is obtained by red-shifting this result ac- 
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FIG. 2: The left panel shows the spectrum of gravitational radiation for a simulation with Lu = 
7vi,/f3, integrated over Ni^. = 32 directions. The error bars result from averaging over eight different 
scenarios. Velocities decrease from top to bottom as Vh = {1, 0.1, 0.01}. The right panel shows 
the parameters A and /*//? as functions of Vb. 
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The two functions and A are displayed in the right panels of Fig. [2] and are approxi- 
mately given by 

0.11 



A 



0.42 + ' 
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(23) 



9 



Kosowsky et al 



This work 



10 



10 



10 



10' 



a 



10 



10 



14 



16 



10 



10 



■20 



1 Mllllll 


1 Mllllll 1 1 1 III 


II 1 1 1 1 null 








/ 


LISA 










--■ BBO 






— , 




/ 

/ 
; 










/ 
/ 

/ 

/ 

/ 

/ 




a 








O 














J 






^1 |»<||||| / 


" 1 


V v\W\ 







10 



10 



10 



10 



12 



10 



10 



16 



10 



10 



10" 



0.01 
f/Hz 



100 



10 



■20 



1 Mllllll 


1 llllllll 1 1 llllllll 1 1 llllllll 






/ 


LISA 




- 




--- BBO 


- 


— ^ 




/ 

/ 

/ 

/ 

; 


- 






/ 

; 








>\ / 












~ 1 / / 


1 1 1 1 im^ 





10 



10' 



0.01 
f/Hz 



100 



FIG. 3: Several spectra of gravitational radiation according to the old and new formulas. The 



parameters are taken from ref. 



and given in table [J with a decreasing from top to bottom. In 



the shaded region, the sensitivity of LISA and BBO is expected to drop considerably. 

Notice that the peak frequency and amplitude agree reasonably well with the results pre- 
sented in ref. {g] (our peak amplitude is about 50% larger for fj, 1, while for small V\, 
both results agree within the statistical errors). In the light of the analysis presented in 



ref. [13|, the dependence of the peak frequency on the wall velocity has the following phys- 
ical interpretation. For small velocities, ff, <^ 1, the phase transition lasts long compared 
to the relevant distance scale that is given by the average bubble size. In this case, the 
GW spectrum inherits the time scale of the source, / ~ /3. If one used in eqs. (l8i)-( |T0l) a 
wall velocity much larger than the speed of light, ff, ^ 1, the phase transition would be 
very short compared to the relevant distance scale and the GWs would inherit the distance 
scale of the source, / ~ This effect leads to a decrease in the peak frequency in the 

transition region where the wall velocity is close to the speed of light. 



IV. CONCLUSIONS 



We reexamined the spectrum of gravitational wave radiation generated by bubble col- 
lisions during a first-order phase transition in the envelope approximation. Using refined 
numerical simulations, our main finding is that the spectrum falls off only as /"^ '^ at high 
frequencies, considerably slower than appreciated in the literature. This behavior is most 
probably related to the many small bubbles nucleated at a later stage of the phase tran- 
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set 


a 


/3/if 


T, 1 GeV 


1 


0.03 


1000 


130 


2 


0.05 


300 


110 


3 


0.07 


100 


85 


4 


0.1 


60 


80 


5 


0.15 


40 


75 


6 


0.2 


30 


70 



TABLE I: Sets of parameters used in Fig. [3j 

sition [3^. This result is especially interesting in the light of recent investigations P, 0] 
that indicate that in the case of a first-order electroweak phase transition (obtained by a 



singlet sector 



3,Q 



or higher dimensional operators 



16 



17 



18| ) the peak frequency of the 



produced radiation is below the best sensitivity range of planed satellite experiments, such 
as LISA and BBO [20I. This effect is shown in Fig. [3] for several typical parameter sets 



for the phase transition in the nMSSM [8|. Notice that the discussion in ref. [8j suggests 
that stronger phase transitions in general lead to smaller peak frequencies due to a decrease 
in the parameters and T*. This amplifies the importance of the high frequency part of 
the gravitational wave spectrum. Notice also that a flatter spectrum simplifies the distinc 



tion from other sources of stochastic gravitational waves, such as turbulence 



or preheating after inflation 25 



26 



27 



28 



21 



23 



24| 



Besides, we found that the peak fre- 
quency slightly depends on the expansion velocity of the bubbles and decreases for higher 
wall velocities. Our quantitative results are summarized by eqs. (fT9l) - (l23l) . 



Finally, we would like to comment on the recent paper [31i, where an analytic approach 
to the GW production by collisions based on stochastic considerations was presented. In 
this approach, assumptions have to be made about the time- dependence of unequal time 
correlations of the velocity field. In their favored model, the authors obtain a scaling as 



for the high frequency part of the spectrum. We suspect that this disparity is due to 



conceptual differences. 

First, notice that the treatment presented here is based on two main ingredients: 
thin wall and the envelope approximations. Even though the stochastic approach in ref. 



'he 



3l| 



does not require the thin wall approximation, the results are also valid in this limit, such that 
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this ap pro ximation cannot be responsible for the different spectra. However, in the approach 



of ref. [31| the colhded and uncoUided regions of the bubbles are treated equally but in a 
stochastic manner. Breaking of the spherical symmetry, necessary for GW production, is 
encoded in assumptions on the velocity correlation functions. This is in contrast to our 
analysis in which the well tested 4] envelope approximation breaks the spherical symmetry 
in a realistic way. 

Furthermore, only the time dependent mean bubble radius enters into the stochastic 
calculation, while in the analysis at hand bubbles with a realistic size distribution are simu- 
lated. The occurrence of many small bubbles probably enhances the high frequency part of 



the spectrum, as already argued in refs. 



Finally, it is interesting to see that even in the stochastic treatment the high frequency 
part of the spectrum can scale as uj~^ (in agreement with our results) if one assumes that the 
source is fully uncorrelated at unequal times. Especially at late stages of the phase transition, 
many small bubbles are generated, start to collide and are absorbed by neighboring bubbles 
at a large rate such that in this regime the assumption of non-correlation might indeed be 
plausible. 
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APPENDIX A: SOME DETAILS ABOUT THE NUMERICAL ALGORITHM 

In this section we will present an efficient way of evaluating the integrals given in eqs. ([8])- 
ffTOj) . It turns out that rotating the vector k parallel to the z-axis is very useful. One 
advantage is that the projection of the transverse-traceless part simplifies to 
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4->4 



FIG. 4: The figure depicts the case of a bubble that intersects with four neighboring bubbles. The 
data structure we use contains the information which bubbles constitute the boundaries of the 
uncollided regions in each segment. This improves the performance greatly especially at the end 
of the phase transition when virtually all bubbles intersect with each other but only a few bubbles 
are relevant in each segment. 



Using cylindrical coordinates and denoting the rotated uncollided region as S'^, this leads to 

dEcw 



Bn,+ {z,t) 



dujdfl 



(1 - Z^' 



AGpl,K^vtu;\\C^\'+\C.\% 



dz e-''"''^^'-'-^'Bn,±{z,t), 



(i0cos(20), Bn-{z,t) 



(A2) 
(A3) 

(A4) 



rf0sin(20). (A5) 



The second advantage of this coordinate system is that the integrals Bn,±{z) do not depend 
on the frequency uj or the wall velocity Vb- Hence, the functions Bn^±{z) need only to be 
determined once for all values of u and Vb what greatly improves the performance of the 
numerical evaluation of the integrals. 

Notice that the function An^ij{u!,t) in eq. (lA4p does not depend on the frequency u and 
the wall velocity Vb separately, but only on the product uvf,. Hence, we choose an exponential 
distribution for the values of the wall velocity Vb and the frequency u to reduce the number 
of necessary evaluations of An^±{uj,t). For each direction, we then evaluate the integrals 
flA2l) - flA5l) and integrate over the different directions. 
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Finally, we would like to comment on the numerical accuracy of our simulation. The 
integration in eq. flA5|) is performed by determining the intersections of the overlapping 
bubbles, while the integration in eq. (]A4I) is done on an adaptive grid such that the relative 
error never exceeds 10^^. This requires for the highest frequencies up to 10^ evaluations of 
the functions Bn,±{z). In order to improve the performance we employ a data structure that 
for a fixed time t and a specific bubble n encodes the regions that are uncollided and do not 
intersect with other bubbles. The z-space is divided into different segments, where in every 
segment the data structure contains the information which bubbles constitutes the borders 
of the intersecting regions. A graphical representation of this data structure is depicted 
in Fig. m This is especially useful at late stages of the phase transition, when almost all 
bubbles intersect with each other, but only a few bubbles are relevant for each segment. 
The dominant error in the final spectrum results from the fact that the time integration 
in eq. (lASO is done on an equidistant grid with size Nt = 512. We checked that the error 
resulting from this choice never exceeds a few percent for the presented results. As a first 
test, we used this algorithm to reproduce the results in the two bubble case and with a cutoff 
function as presented in ref. jsj. 

A final source of uncertainty is the finite size of the volume. As argued in the main 
text, the completion of the phase transition takes approximately the time r ~ 5//5 and the 
the size of the volume under consideration should be chosen accordingly. It turns out that 
for the portion of the spectrum under consideration, a size of Lu = Ivb/ (3 is sufficient to 
suppress boundary effects, as long as the first bubble nucleates not too far from the center 
of the volume which we assure in our simulations. 
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